MA5701 Optimización no lineal: Tarea 3

Fecha de entrega: 4 de Julio a las 23:59 hrs.

Profesor: Jorge G. Amaya A.

Auxiliar: Aldo Gutiérrez Concha.

Ayudantes: Carolina Chiu y Mariano Vazquez.

Autor: Felipe Urrutia Vargas

Repositorio: github.com/furrutiav/nonlinear_optimization_zoutendijk

Tabla de contenidos

Método de direcciones admisibles (Zoutendijk)

En este trabajo nos enfocaremos en resolver el problema de optimización no lineal: $$(P) \quad \begin{cases} \min &f (x) \\ s.a. &Ax \leq b \\ &Ex = e \end{cases}$$ en su generalidad.

Para esto, se realizara una implementacion del metodo Método de direcciones admisibles o metodo de Zoutendijk. Algoritmo que sigue los siguientes pasos:

Pseudo-algoritmo

(0) Sean $\epsilon>0$, $k = 0$, $x_0 \in \mathbb{R}^n$ tal que $Ax_0 \leq b$, $Ex_0=e$

(1) Sea la descomposición $A = [A_1,A_2]^T, b =(b_1, b_2)^T$ tal que $A_1x_k = b_1, A_2x_k < b_2$.

(2) Resolver el problema lineal

$$(D_k) \quad \begin{cases} \min &\nabla f (x_k)^T d\\ s.a. &A_1d \leq 0 \\ &Ed = 0\\ &−1 \leq d_j \leq 1, \quad j = 1, . . . , n \end{cases}$$

y sea $d_k$ solución de $(D_k)$.

Si $|| \nabla f (x_k)^T d_k || < \epsilon$, entonces parar.

Si no, ir a (3).

(3) Determinar el paso, resolviendo aproximadamente el problema de minimización unidimensional

$$(L) \quad \begin{cases} \min &f (x_k + \lambda d_k)\\ s.a. &\lambda \in [0, \bar{\lambda_k}] \end{cases}$$

mediante el método de Armijo.

Se usa

$$ \bar{\lambda_k} = \min \left\lbrace \frac{(b_2 − A_2 x_k)_i}{(A_2 d_k)_i} : (A_2 d_k)_i > 0 \right\rbrace$$

y se considera $\bar{\lambda_k} = +\infty$ cuando $(A_2 d_k)_i \leq 0, \forall i$.

Sea $\lambda_k$ la solución del subproblema $(L)$. Hacer:

$$x_k+1 = x_k + \lambda_k d_k,$$$$k ← k + 1,$$

e ir a $(1)$.

Implementación

Contenidos

Esta seccion se divide en los siguientes capitulos:

Los capitulos Paso i corresponden a implementar el i-esimo paso del Pseudo-algoritmo.

Mientras que el capitulo Clase Zoutendijk corresponde al desarrollo de una clase (paradigma de programacion de python) que permite ingresar y resolver cualquier problema de tipo $(P)$.

Librerias

La siguiente celda contiene todas las librerias utilizadas en este trabajo.

Paso 0

Los parametros del paso 0 son $\epsilon>0$, $k = 0$, $x_0 \in \mathbb{R}^n$ tal que $Ax_0 \leq b$, $Ex_0=e$. Notar que $\epsilon$ y $k$ son arbitrarios. Por otro lado, el punto inicial $x_0$ debe de satisfacer las restricciones lineales $Ax\leq b$ y $Ex=e$. Para verificar la factibilidad del punto inicial, se utiliza la siguiente función:

Paso 1

Dado un punto factible $x_k$, el objetivo del paso 1 es descomponer las restricciones de desigualdad por componentes activas e inactivas. Es decir, una descomposición de $A = [A_1,A_2]^T$ y $b =(b_1, b_2)^T$ tal que $A_1x_k = b_1, A_2x_k < b_2$. Para resolver esto, se utiliza la siguiente función:

Paso 2

Dado un punto factible $x_k$ tal que $A_1 x_k = b_1$, el objetivo del paso 2 es resolver el siguiente problema lineal:

$$(D_k) \quad \begin{cases} \min &\nabla f (x_k)^T d\\ s.a. &A_1d \leq 0 \\ &Ed = 0\\ &−1 \leq d_j \leq 1, \quad j = 1, . . . , n \end{cases}$$

Para determinar el gradiente de la funcion en el punto $x_k$, utilizamos la libreria $\texttt{numdifftools}$ con la clase $\texttt{Gradient}$. Para leer la documentacion, revisar aqui.

Para resolver $(D_k)$ en python, consideramos la libreria $\texttt{scipy}$ con el modulo $\texttt{optimize}$. Entra las funciones de optimizacion disponibles, utilizaremos la funcion $\texttt{linprog}$ que se encarga de resolver este tipo de problema de programacion lineal. Para mas detalle, revisar la documentacion disponible aqui. Por defecto, utilizaremos el método simplex para resolver cada problema $(D_k)$.

Llamaremos por $d_k$ a la solución de $(D_k)$. Para la siguiente proposicion $| \nabla f (x_k)^T d_k | < \epsilon$ (criterio de parada), avisaremos con una variable $\texttt{stop}$ su valor de verdad (verdadero o falso).

En la celda que sigue, se encuentra una funcion que resuelve $(D_k)$ y entrega tanto la direccion $d_k$ encontrada como el valor de verdad, del criterio de parada, en una variable $\texttt{stop}$:

Paso 3

Dada un direccion de descenso $d_k$ en un punto factible $x_k$, el objetivo del paso 3 es determinar el paso $\lambda_k$. Para esto, se propone resolver aproximadamente el siguiente problema de minimización unidimensional:

$$(L) \quad \begin{cases} \min &f (x_k + \lambda d_k)\\ s.a. &\lambda \in [0, \bar{\lambda_k}] \end{cases}$$

Para determinar $\bar{\lambda_k}$, se usa

$$ \bar{\lambda_k} = \min \left\lbrace \frac{(b_2 − A_2 x_k)_i}{(A_2 d_k)_i} : (A_2 d_k)_i > 0 \right\rbrace$$

y, en caso de que el conjunto sea vacio, se considera $\bar{\lambda_k} = +\infty$. En la implementacion, $+\infty = 1$.

Parte de este trabajo es resolver $(L)$ mediante el método de Armijo. Este metodo consiste en lo siguiente: Se define $\varphi(t) = f (x_k + t d_k)$ y se elige $\sigma \in (0,1)$ de tal forma que la recta de pendiente $\sigma \varphi'(0)$ se utiliza para acotar superiormente a la funcion en una vecindad del cero; y sobre este dominio encontrar el paso mas largo posible. En efecto, encontrar

$$t^* = \max \{t : \varphi(t) \leq \varphi(0) + \sigma \varphi'(0)t\}$$

Para resolver este otro problema, consideramos la siguiente aproximacion:

$$t^* = \max \{t : \exists p\in \mathbb{N}, t=ph, \varphi(t) \leq \varphi(0) + \sigma \varphi'(0)t\}$$

donde $h>0$ es dado por el usuario. En este trabajo, se definira $h = \frac{\bar{\lambda_k}}{N}$ con $N>0$ grande y dado por el usuario.

Para determinar la derivada de la funcion $\varphi$ en el cero, utilizamos la libreria $\texttt{numdifftools}$ con la clase $\texttt{Gradient}$. Para leer la documentacion, revisar aqui.

Sea $\lambda_k$ la solución del subproblema $(L)$. Con dicha solucion, hay que generar la siguiente iteracion dado el paso y la direccion:

$$x_{k+1} = x_k + \lambda_k d_k,$$

Finalmente, en la celda que sigue se encuentra la funcion que resuelve el problema $(L)$, y retorna el nuevo punto $(x_{k+1})$. Adicionalmente, entrega al valor de verdad de la siguiente proposicion $\lambda_k< \epsilon$. Dicha proposicion, es otro criterio de parada que consideramos produente para el algoritmo.

Clase Zoutendijk

En esta seccion se realiza un ensamble de los pasos ya definidos. Adicionalmente, se diseña un metodo denominado $\texttt{solve}$ que ejecuta este ensamble dado ciertos parametros iniciales. Dicho metodo, es parte de la clase principal de esta seccion que denominaremos como $\texttt{Zoutendijk}$.

La idea es que esta clase, logre resolver en su generalidad cualquier problema de tipo $(P)$. Para esto, se añade un lector con expresiones regulares para leer las restricciones. De tal forma que dado un string logre generar las metrices de desigualdades e igualdades automaticamente. Asi mismo, la funciones objetivo es escrita en formato $\texttt{python}$. Por ejemplo, la funcion $x^2$ se introduce por $\texttt{x**2}$.

Ahora bien, consideramos los siguientes criterios de parada para el metodo que resuelve los problemas de tipo $(P)$. Estos criterios son:

Cada uno logra evitar posibles situaciones donde realizar otra iteraciones es marginal al resultado final. Se añade el numero maximo de iteraciones para evitar casos patologicos como por ejemplo entrar en un loop que generen iteraciones ciclicas. Por otro lado, detenerse si la siguiente iteracion es infactible, caso que puede ocurrir por errores numericos. La restriccion de $||\nabla f(x_k)^T d_k|| < \epsilon$ es la que viene por defecto con el metodo de $\texttt{Zoutendijk}$. Sin embargo, se considero tambien $\lambda_k < \epsilon$, que se traduce en un paso pequeño.

En lo que sigue se encuentra la clase ya descrita. Por simplicidad, pasaremos directamente a los resultados:

Resolución de problemas

En esta seccion se realizara la resoluciones de dos problemas con el metodo de Zoutendijk.

Problema 1

$$(P_1) \quad \begin{cases} \min &8(x_1-6)^2+(x_2-2)^4 \\ s.a. &-x_1+2x_2 \leq 4 \\ &3x_1+2x_2 \leq 12 \\ &x_i \geq 0, \quad i=1,2 \\ \end{cases}$$

Para comparar nuestros resultados con el optimo real del problema $(P_1)$ consideramos el motor https://www.wolframalpha.com. La instancia de este problema se puede encontrar aqui. Que entrega la siguiente solución:

x1 x2 f(x1, x2)
3.858319799174562844633020514 0.2125203012381557330504692291 46.90291234143656803690204873

Inicializar datos del problema

Este es un problema de dos variables con cuatro restricciones. Todas ellas de desigualdad.

Resolver problema

Inicializamos con un punto incial, epsilon pequeño, sigma igual a $0.2$ y un $N=100$. Donde $N$ es tal que $h_k = \frac{\bar{\lambda}_k}{N}$ es el $h$ del paso de Armijo. Es decir, cuando $N$ es mas grande entonces $h_k$ es mas pequeño. Ademas, consideramos un maximo de $100$ iteraciones.

Del reporte final obtenemos que el metodo tuvo que detenerse por exceder el maximo de iteraciones.

Reportar resultados

Error c/r solución

En la siguiente parte calculamos el error entre el punto encontrado y la solucion, y el error entre los valores de la funcion objetivo en la ultima iteracion y el valor del problema.

Podemos observar que los errores son altos, verificando que el punto encontrado no es la solucion real

GridSearch

Dado que los resultados obtenidos no fueron los mejores, realizaremos un analisis de los hiperparametros del metodo para entender cual es una configuracion adecuada para resolver este problema. Para esto, solamente probaremos distintos valores de N y sigma. Para el primero se estudiaran valores de la forma $10^k$, con $k=1,...,5$. Y para el segundo, $21$ valores en $[1e6, 0.9]$. Con esto obtenemos una grilla de 105 pares de valores diferentes. Esto se traduce en 105 instancias diferentes de hiperparametros para resolver el problema 1.

Graficaremos la grilla y pintaremos las instancias segun el criterio de parada que utilizo el metodo para detenerse.

Podemos observar que para valores de sigma menores que 0.4 la mayoria de las instancias se detienen por superar el maximo numero de iteraciones. Por otro lado, existe un porcentaje no menor de instancias que se detienen por obtener iteraciones infactibles. Para estos casos, no existe un patron claro. Por otro lado, el resto de instancias se detiene por un paso pequeño.

Graficamos solo las instancias que son factibles. Y se colorea la grilla segun el error obtenido comparando la funcion objetivo en la ultima iteracion con el valor real del problema 1.

Se observa que en las mismas instancias en que el metodo se detuvo por superar el maximo de iteraciones, entonces tambien tiene un error alto.

De lo observado anteriormente, solo visualizaremos aquellas instancias que se detuvieron ya que el paso era muy pequeño. Esta vez, el eje x es el valor del problema, el eje y es el valor de sigma y se colorea cada instancia segun el N.

Del grafico anterior, para casi cualquier valor de sigma, se observa que si se elige un N grande entonces, el metodo logra obtener soluciones mas cercanas a la solucion real.

Mismo grafico anterior, pero esta vez los colores son en funcion al tiempo de ejecucion.

Del grafico anterior no es claro si los parametros afectan mas que otros para obtener mejores soluciones en menor tiempo.

De la busqueda de hiperparametro realizada anteriormente, se eligen valores de N=$10^5$ y sigma igual a $0.495$. Para lo cual, se obtienen los siguientes resultados:

En solo 11 iteraciones el metodo logra llegar al optimo.

Problema 2

$$(P_2) \quad \begin{cases} \min &x_1^4-2x_2^2+10x_1x_2^2+x_4^2 \\ s.a. &x_1+x_2-x_3 = 1 \\ &x_1+x_4 = 4 \\ &x_1-x_2=0 \\ &x_i \geq 0, \quad i=1,2,3,4 \\ \end{cases}$$

Para comparar nuestros resultados con el optimo real del problema $(P_2)$ consideramos el motor https://www.wolframalpha.com. Para esto, utilizamos una version simplificada, dado que dicho problema se puede transformar en uno unidimensional. Esto es:

$$(P_2) \quad \begin{cases} \min &x_1^4+10x_1^3-x_1^2-8x_1+16\\ s.a. &x_1\geq 1/2\\ &x_1\leq 4\\ &x_2=x_1\\ &x_3=2x_1-1\\ &x_4=4-x_1 \end{cases}$$

La instancia de este problema se puede encontrar aqui. Que entrega la siguiente solución para $x_1$

x1 f(x1)
0.5311288741492748261833066 13.04675363940572740382986

Inicializar datos del problema

Este es un problema de 4 variables, con 7 restricciones. 3 de igualdad y 4 de desigualdad.

Resolver problema

Inicializamos con un punto incial, epsilon pequeño, sigma igual a $0.2$ y un $N=100$. Donde $N$ es tal que $h_k = \frac{\bar{\lambda}_k}{N}$ es el $h$ del paso de Armijo. Es decir, cuando $N$ es mas grande entonces $h_k$ es mas pequeño. Ademas, consideramos un maximo de $100$ iteraciones.

Del reporte final obtenemos que el metodo tuvo que detenerse por exceder el maximo de iteraciones.

Reportar resultados

Error c/r solución

En la siguiente parte calculamos el error entre el punto encontrado y la solucion, y el error entre los valores de la funcion objetivo en la ultima iteracion y el valor del problema.

Podemos observar que los errores son altos, verificando que el punto encontrado no es la solucion real

GridSearch

Al igual que los resultados obtenidos preliminarmente para el problema 1, los resultados obtenidos para el problema 2 no fueron los mejores. Por esta razon, realizaremos un analisis de los hiperparametros del metodo para entender cual es una configuracion adecuada para resolver este problema. Para esto, solamente probaremos distintos valores de N y sigma. Para el primero se estudiaran valores de la forma $10^k$, con $k=1,...,5$. Y para el segundo, $21$ valores en $[1e6, 0.9]$. Con esto obtenemos una grilla de 105 pares de valores diferentes. Esto se traduce en 105 instancias diferentes de hiperparametros para resolver el problema 2.

Graficaremos la grilla y pintaremos las instancias segun el criterio de parada que utilizo el metodo para detenerse.

Podemos observar que para valores de sigma menores que 0.4 la mayoria de las instancias se detienen por superar el maximo numero de iteraciones. Por otro lado, existe solo una instancia que se detienen por el criterio del gradiente producto la direccion es pequeño. Por otro lado, el resto de instancias se detiene por un paso pequeño.

Graficamos todas las instancias que son factibles y se colorea la grilla segun el error obtenido comparando la funcion objetivo en la ultima iteracion con el valor real del problema 2.

Se observa que en las mismas instancias en que el metodo se detuvo por superar el maximo de iteraciones, entonces tambien tiene un error alto.

De lo observado anteriormente, solo visualizaremos aquellas instancias que se detuvieron ya que el paso era muy pequeño. Esta vez, el eje x es el valor del problema, el eje y es el valor de sigma y se colorea cada instancia segun el N.

Del grafico anterior, para casi cualquier valor de sigma, se observa que si se elige un N grande entonces, el metodo logra obtener soluciones mas cercanas a la solucion real.

Mismo grafico anterior, pero esta vez los colores son en funcion al tiempo de ejecucion.

Del grafico anterior no es claro si los parametros afectan mas que otros para obtener mejores soluciones en menor tiempo.

De la busqueda de hiperparametro realizada anteriormente, se eligen valores de N=$10^4$ y sigma igual a $0.495$. Para lo cual, se obtienen los siguientes resultados:

En solo 6 iteraciones el metodo logra llegar al optimo.

Conclusión

En este trabajo nos hemos propuesto resolver en general problemas de optimización de tipo $(P)$ con restricciones lineales, tanto de desigualdad como de igualdad. Para resolver este tipo de problemas definimos el pseudoalgoritmo de las direcciones admisibles (o método de Zoutendijk).

Para implementar este método, las funciones que resuelven cada uno de los pasos del algoritmo se definen de forma independiente. Posteriormente, los pasos se ensamblan en una única función como método de una clase de Python, llamada $\texttt{Zoutendijk}$. Además de resolver problemas de tipo $(P)$, es capaz de leer automáticamente la función objetivo y las restricciones del problema a partir de un texto de estilo Python.

Se propusieron dos problemas de optimización para resolver con Zoutendijk. Cuyas soluciones se obtuvieron satisfactoriamente a partir de una búsqueda de hiperparámetros. Para cada problema se estudió el comportamiento del método para diferentes valores de N y sigma. Parámetros que son clave en el Paso de Armijo. El primero modula lo fino que es el dominio de búsqueda y el segundo caracteriza el dominio a partir de un recta. Para encontrar el mejor par de hiperparámetros para cada problema se utilizó la técnica GridSearch.

Del análisis de los parámetros obtuvimos algunas conclusiones interesantes. Entre estos estan: (1) Valores pequeños de sigma (menos de 0,4) hacen que el método utilice demasiadas iteraciones para alcanzar incluso soluciones no optimales, (2) Para valores mayores de sigma (mayores de 0. 4) el método consigue utilizar menos de 100 iteraciones y detenerse antes de tiempo debido a un paso pequeño, (3) Utilizar valores grandes de N (igual a 10^5 para el problema 1 y =10^4 para el problema 2) permite obtener soluciones más cercanas a la solución real, y (4) No se obtuvieron patrones claros para el tiempo de ejecución del código, sin embargo resolver cada problema no toma más de 1 segundo.

Finalmente, todo este trabajo está disponible en un repositorio de github. Véase el repositorio en github.com/furrutiav/nonlinear_optimization_zoutendijk. Además, quedan experimentos futuros para mejorar el rendimiento del método. Esto consiste en mejorar el método de Armijo con inicializaciones aleatorias de las variables y muestreo aleatorio del dominio de búsqueda. Por otro lado, queda por revisar el rendimiento del método para otros tipos de funciones objetivo que no sean polinómicas y regiones de decisión más intrincadas. Entre ellas se encuentra la conocida función Rosenbrock.